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ABSTRACT 

Recent observational studies of type la supernovae (SNela) suggest correlations between the peak 
brightness of an event and the age of the progenitor stellar population. This trend likely follows from 
properties of the progenitor white dwarf (WD), such as central density, that follow from properties 
of the host stellar population. We present a statistically well-controlled, systematic study utilizing 
a suite of multi-dimensional SNela simulations investigating the influence of central density of the 
progenitor WD on the production of Fe-group material, particularly radioactive 56 Ni, which powers 
the light curve. We find that on average, as the progenitor's central density increases, production of 
Fe-group material does not change but production of 56 Ni decreases. We attribute this result to a 
higher rate of neutronization at higher density. The central density of the progenitor is determined 
by the mass of the WD and the cooling time prior to the onset of mass transfer from the companion, 
as well as the subsequent accretion heating and neutrino losses. The dependence of this density 
on cooling time, combined with the result of our central density study, offers an explanation for the 
observed age-luminosity correlation: a longer cooling time raises the central density at ignition thereby 
producing less 56 Ni and thus a dimmer event. While our ensemble of results demonstrates a significant 
trend, we find considerable variation between realizations, indicating the necessity for averaging over 
an ensemble of simulations to demonstrate a statistically significant result. 

Subject headings: hydrodynamics — nuclear reactions, nucleosynthesis, abundances — supernovae: 
general — white dwarfs 



1. INTRODUCTION 

Observations targeting the environment of type la 
supernovae (SNela) have exposed open questions con- 
cerning the dependence of bot h their rates and aver - 
age brightness on environment. iMannucci et al.l (2006) 
show that the dependence of the SNIa rate on delay 
time (elapsed time between star formation and the su- 
pernova event) is best fit by a bimodal delay time dis- 
tribution (DTD) with a prompt component less than 
1 Gyr after star formation and a tardy component sev- 
eral Gyr later. Although the clar ity of this effect is 
clouded by galaxy sampling (JFihppcnko 200 9J), the ba- 
sic re s ult is borne out even w ithin galaxies (jRaskin et al.l 
I2009T ). iGallagher et al.l (J2008I ) measure a correlation be- 
tween the brightness of a SNIa and its delay time, which 
they state is consistent with either a bimodal or a contin- 
uous DTD Other recen t studies bvlHowell et al.l ([2009), 
iNeill et al.l (|2009fl . and iBrandt et alj (|2010T> also find a 
correlation between the delay time and brightness of a 
S NIa. 

Phillips (1993) identified a linear relationship between 
the maximum B-band magnitude of a light curve and 
its rate of decline. This "brighter equals broader" re- 
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lationship has been extended to additional bands with 
templates from nearby events, allowing SNela to be cali- 
brated as an extension of the astronomical distance lad- 
der fsee lJha et al.ll2007l for a description of one method). 
The brightness of a SNIa is determined principally by the 
radio active decay of 56 Ni synthesized during the explo- 
sion dTruran et alJll967tlColgate fc McKed ll969l: lArnettl 
1982: Pinto fc Eastmanll2000f) . 

A widely-accepted proposal to explain many, if not 
most, events is the thermonuclear disruption of a white 
dwarf (WD) in a mass-transferring binar y system (for re- 
views from various perspectives see: [Branch et al.l|199E 



Filippcnk d 119971: IHillebrandt k Niemeverl 120001: ILiviol 
2000; Ropke 2006). In this paradigm, a longer delay 
time suggests the possibility of a longer elapsed time be- 
tween the formation of the WD and the onset of accre- 
tion. During this period, denoted here as the WD cooling 
time (tcooi), the WD is isolated from any significant heat 
input and decreases in temperature. A longer £ coo i re- 
sults in a higher central density when the cor e reaches 
the ignition temperature (Lesaffre et al. 2006), due to 
the lower entropy at the onset of accretion. Thus, a cor- 
relation between central density and the peak brightness 
of an event suggests a correlation between delay time 
and the brightness of an event. While previous work in- 
dicated a correlation between central density and peak 
brightness, none has averaged o yer a statistically sig- 
nificant ensemble of realizations dBrachwitz et al.|[2000t 
IRopke et al.|[200l IHoflich et al.H201Q[) . Therefore, we in- 
vestigate, for the first time, a statistically significant cor- 
relation between progenitor central density and average 
peak brightness of SNela. 



The surrounding stellar population, the metallicity 
and mass of the progenitor, the thermodynamic state 
of the progenitor, the cooling and accretion history of 
the progenitor, and other parameters are known to af- 
fect the light curves of SNela; the role, and even pri- 
macy, of these va rious parameters i s the subject of on- 
going study (e.g..lRopke et all 120061: iHoflich et alll2010t 
Uackson et al.H2010D . Additionally, many of these effects 
may be interconnected in complex ways (jLesaffre et al.l 
2006) . In this study, we isolate the direct effect of vary- 
ing the progenitor central density on the production of 
56 Ni. To first order, this yield controls the brightness of 
an event; second-order effects on the light curve are left 
for future study. 

2. METHOD 



Once a WD forms in a binary, it is initially isolated and 
slowly cools in a single-degenerate scenario. Eventually, 
mass-transfer begins to carry light elements from the en- 
velope of the companion to the surface of the WD. If 
the accretion rate exceeds a thresh old, the infalling ma - 
terial experiences steady burning ()Nomoto et al.|[2007f ). 
and eventually the WD gains enough mass to compress 
and heat the core. Once the temperature rises enough 
to initiate carbon reactions, the core begins to convect 
(or "simmer" ) . Our progenitor models parameterize the 
WD at the end of this simmering phase, just prior to the 
birth of the flame that eventually will disrupt the entire 
WD in a SNIa. 

We constructed a series of five parameterized, hy- 
drostatic progenitor models that account for simmer- 
ing in which we vary the central density (p c ). The 
outer regions are isothermal, although some te mper- 
ature structure is expected ()Kuhlen et al.l [2006) , and 
the cores are ise ntropic due to convection and have a 
lower C/O ratio dStraniero et all 120031: iPiro fe Bildstenl 
[200l iChamulaket all 120081 : IPiro fe Changl 12008) . 
(jJackson et al.ll2010f ) explored these effects and we chose 
the core composition as 40% 12 C, 57% 16 0, and 3% 
22 Ne and the outer layer as 50% 12 C, 48% 16 0, and 
2% 22 Ne. For our p c we chose 1 — 5 x 10 9 g/cm 3 in 
steps of 10 9 g/cm 3 . The central temperature must be 
in the range of carbo n ignition, which is approximately 
7 - 8 x 10 8 K (e.g., iKuhlenet allllZOOll : we selected 
7 x 10 8 K. Based on prior research, we chose other 
model parameters to produce exp ected amounts of Fe- 
group elements in the explosion (Tow nsley et al„ 2009; 
Uackson et al.l[2010h . The values were kept constant in 
all simulations in order to isolate the central density ef- 
fects. 

With these five progenitor models, we utilize the st atis- 
tical framework presented in iTownslev et al.l (|2009l ) for 
a controlled study of the effect of varying the central 
density. For each progenitor, we created thirty realiza- 
tions seeded by a random number used to generate a 
unique set of spherical harmonics with power in modes 
12 < £ < 16. The spectra are used as initial pertur- 
bations to the spherical flame surface around the center 
of the progenitor star. Each progenitor uses the same 
seed values, resulting in the same thirty perturbations. 
This choice allows us to check for systematic biases in 
the realizations across different progenitors. 

We performed a suite of 150 two-dimensional, ax- 
isymmetric simulations of the deflagration-to-detonation 



transition (DDT) model of SNela with a customized 
version of FLASH, a compressible, Eulerian, adaptive- 
mesh, hydrodynamics code. The modifications to this 
code are (1) the burning model, (2) the flame speed 
computations, (3) the mesh refinement criteria, (4) the 
DDT criterion. Our simulation me thods are described 
in detail in previous public ations (|C alder et al.l 120071 : 

md continue to b e im- 
proved (jJackson_et_aL| [20101; iTowHiixiEin HOIS) • We 



ITownslev et~al] 12 007. 



pupiic atu 
, [2001) ai 



should note that the DDT criterion is based on a char- 
acteristic density at which we ignite detonations, which 
we select to be 10 71 g/cm 3 Fo llowin g the procedures 
descri bed in I C alder et al.l (|2007l ) and ISeitenzahl et al.l 
(2009) for calculating the neutronization rate in material 
in nuclear statistical equilibrium (NSE), we utilize weak 
rates from iFuller et all (119851 ) . IQda et al.l (|1994| ). and 
Langanke & Martincz-Pincdo (2001), with newer rates 
superseding earlier ones. The reaction networks for cal- 
culating the energetics and time scales of the deflagra- 
tion and detonation phases included the same 200 nu- 
clides, and the NSE calculation of tables included 443 
nuclides. The supernova simulations use a three-stage 
burning model in which the timescale to burn to NSE 
is calibrated to reproduce the correct yield of Fe-group 
elements. The weak reaction rate is negligible at a den- 
sity where Si-burning to 56 Ni is incomplete; therefore, 
we estimate the 56 Ni yield from the total NSE yield 
and its associated electron fraction. We performed nu- 
clear postprocessing of Lagrangian tracer particles us- 
ing a network of 200 nuclides on a subset of simula- 
tions. Comparison demonstrates that our burning model 
matches the energetics a nd final composition of impor- 
tant nuclei such as 56 Ni (ITownslev et al.H2010D (see als o 
studies bv lTravaglio et al.ll2004l: ISeitenzahl et al.ll2010f ). 
We utilize the adaptive-mesh capability, with a highest 
resol ution of 4 km, which demonstrates a converged re- 
sult (|Townslev et al.ll2009t ). 

3. RESULTS AND DISCUSSION 

The 56 Ni yield is determined partly by neutronization 
occurring during the thermonuclear burning. Neutron- 
ization pushes the nucleosynthetic yield away from bal- 
anced nuclei such as 56 Ni to more neutron-rich, stable 
isotopes like 58 Ni. Thus the amount of neutronization 
influences the brightness of an event and, all else be- 
ing constant, more neutronization results in a dimmer 
event. The degree of neutronization depends on the den- 
sity and temperature evolution of burned material. Gen- 
erally, thermonuclear burning occurring at higher densi- 
ties will neutronize faster. In an explosive event like a 
supernova, the longer material remains i n NSE at high 
densities, the more neutronization o ccurs (Nomot o et al.1 
[J98i IKhokhTovl fl99ll ICalder et al.ll2007h . Accordingly, 
for SNela, both the central density and the duration 
of the deflagration phase influence the brightness of an 
event. 

Figure [T] presents the mass fraction of NSE material 
that is 56 Ni as a function of deflagration duration, with 
points colored to indicate p c . The duration of the defla- 
gration phase is the time elapsed between the formation 
of a flame front and the ignition of the first detonation 
point. We consider this elapsed time because there is 
little contribution to neutronization after the first DDT. 
The mass of NSE material produced increases during the 
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Fig. 1. — Plot of the 56 Ni-to-NSE mass ratio vs. duration of 
deflagration. Small crosses are single simulations and are colored 
by p c . Large black crosses show the average values for a given p c 
with error bars showing the standard deviation and the standard 
error of the mean. 



course of the SNIa and eventually plateaus; we find the 
point where the NSE yield changes by less than 0.01% 
over 0.01 s and use that mass as the final yield of the 
SNIa. The results have considerable scatter but show two 
trends. First, at a given p c , simulations with longer defla- 
gration periods tend to have a greater degree of neutron- 
ization. Next, simulations from progenitors with higher 
central densities tend to have a greater degree of neutron- 
ization despite having shorter deflagration periods. This 
result shows that the increased rate of neutronization at 
higher densities, which can be seen by the steeper slopes 
of the higher density trend lines, more than compensates 
for the decrease in time for neutronization to occur. Ac- 
cordingly, our results qualitatively agree with previous 
work (e.g., Ilwamoto et al.l ll999t iBrachwitz et all 120001: 
IHofhch et al.ll2010fl . We~also note that the yield of NSE 
material is, with in the error of the s lope, independent 
of p c . As seen in lWooslevet all (2007), if the NSE yield 
remains constant but the amount of 56 Ni varies, then the 
results should lie along the Phillips relation. 

The greater degree of neutronization seen in Figure [1] 
leads directly to lower 56 Ni yields with increasing p c . 
This result can be seen in Figure [3J which presents the 
56 Ni yield of each simulation plotted against p c with color 
used to classify the simulations by realization. The scat- 
ter among different realizations at the same p c is greater 
than the variation across the p c range, indicating the 
need for analysis of an ensemble of realizations. This 
scatter can result in a single realization showing a trend 
unlike the statistical trend; for example, by considering 
only realization 2 and p c of 2 x 10 9 and 3 x 10 9 g/cm 3 
(Figure O green curve), we would conclude that increas- 
ing p c causes an increase in 56 Ni production, instead of a 
decrease as seen in the overall ensemble. The scatter fol- 
lows from a strong dependence on the morphology of the 
flame surface during the early deflagration, which varies 
the duration of the deflagration and the production of 
56 Ni. Changing p c can cause a local change in the plume 
dynamics that overrides the general trend. By fitting the 
averages we find the relation 
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Fig. 2. — Plot of mass of 56 Ni produced vs. p c for five different 
realizations (colored curves), demonstrating the variety of trends 
seen for a single realization. These include non-monotonic trends, 
which could suggest an increase of 56 Ni with increasing p c instead 
of the decrease seen in the ensemble. In black are the average values 
for each density, along with the standard deviation, the standard 
error of the mean, and a regression fit to the average values. 



where 



Mr. 



,4 = -0.047 ±0.003 

10 9 g/ cm 
B = 0.959 ±0.009 M . 



AL 



5B Ni 



A Pc + B, 



(1) 



We use data from lLesaffre et ail (|2006T ) to correlate p c 
to icooi, one component of the delay time, specifically 
the results for a WD with a pre-accretion mass of 1 Mq . 
Thus if we imagine a collection of stars forming with 
the same zero-age main-sequence mass but with differ- 
ent companions and binary separations, the delay time 
would be dominated by i coo i. There are large uncertain- 
ties in the relationship between p c and t C ooi; accordingly, 
we neglect error derived from this la rge uncertainty in our 
analysis. We note that the work of iLesaffre et al.l (2006) 
suggests that a WD with a central density of 10 9 g/cm 3 
will not ignite; further accretion is necessa ry to reach ig- 
nition conditions. Therefore we cannot use ILesaffre et all 
(2006) to compute a t coo i for our lowest-p c simulations, 
and so we omit such simulations. 

We convert our 56 Ni mass to stretch (s) (jHowell et al.l 
2009), the magnitude difference in the B band be- 
tween maximum l ight and 15 days aft er maximum 
light (Am 15 (B)) (jGoldhaber et al.l 1200 ID . and abso- 
lute mag nitude in the V band at maximum light 
(My) (Phillips et al.lll999l ) in order to compare with ob- 
servational findings, shown in Figure [31 A detailed com- 
parison of our models to these observables requires ra- 
diative transport calculations of synthetic spectra and 
light curves, but our results drawn from the 56 Ni mass 
are sufficient f o r the basic properties we consider here. 
iWooslev et al.l (|2007f ) showed that the brightness of an 
event depends on the mass of 56 Ni when it is distributed 
through a large fr action of the star as in our simula- 
tions. Additionally, iStritzinger et al.l (J2006D showed that 
late-time nebular spectroscopy finds 56 Ni yields consis- 
tent with those found from peak luminosity using the 
inverse of the relations we use. 



-19.5 




1 i ' ' i 


f -19.4 










- 


-19.3 










1.3 




1.25 






" 


5 1 - 2 










- 


1.15 














1.1 


^ 


2.0e9g/crru i — i — i 
3.0e9g/crrr i — i — i 
4.0e9g/cm, i — i — i 
5.0e9 g/crri i — i — i 




0.9 


T 




£ 0.8 

ID 

I 0.7 














0.6 




■ i ■ ■ i 



10" 



10° 



white dwarf cooling time (yr) 



Fig. 3. — Plot of brightness vs. t coo ] in terms of My (upper 
panel), s (central panel), and Amis(B) (lower panel). Plotted are 
the average values for each p c along with the standard deviation, 
the standard error of the mean, and a regression fit to the average 
values. 



TABLE 1 
Best-fit parameters for brightness-age relations. 



Ai 



M v 



AmisCg) 



0.078 ± 0.006 
-0.064 ± 0.005 
0.120 ±0.009 



-20.05 ±0.05 

1.74 ±0.04 

-0.26 ±0.08 



We find that the best-fit relations follow the form 



q = a q log! 



yr 



■0 q 



(2) 



where q is one of Ami^(B), s, or My. The values for a q 
and /3q are shown in Table [T] 

To highlight the comparison with observations, in Fig- 
ure 2] we plot an expansion of the center pane l from Fig- 
ure[3]with the binned results from Figure 5 of Ncill et al. 
(2009). While an absolute comparison is not possible, the 
similarity of the overall trend indicates that variation of 
p c is an important contribution to the observed depen- 
dence. Our choice of initial conditions and DDT den- 
sity results in an effective calibration that yields higher 
than expected 56 Ni masses. Accordingly, our results are 
systematically too bright, giving abnormally high val- 
ues of stretch. Future studies will correct for this ef- 
fect. A more sub tle point comes in the usage of "age": 
iNeill et all (|2009f ) measures the luminosity- weighted stel- 
lar age, while for the theory we have simply used t coo i 
directly, which for late times is the dominant portion of 
the time elapsed since star formation. Such offsets, ei- 
ther vertical or horizontal, are less important than the 
overall trend and the range that can be attributed to 
variation of p c . For comparison, the black line in Fig- 
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Fig. 4. — Plot of stretch vs. age. In red are the points from this 
study, based on variations in p c , along with the standard deviation, 
the standard error of the mean, and a best-fit trend line following 
the form of Equation J2Jj. In blue are th e binned and averaged 
points from Figure 5 of Ncill ct al. (2009), along with a best-fit 
trend line following the form of Equation J2}. The vertical grey 
lines mark the cuts between bins in the Ncill ct al. (2009) analysis. 
The trend line for the Neill ct al. (2009) data is shifted upward for 
comparison to our results (black line). The overall offset to larger 
stretch in the simulations is due to the choice of DDT density. 
The approximate agreement of the overall trend indicates that the 
variation of p c is an important contributor to the observed trend, 
but that other factors are also important. 



ure[H shifts the best-fit line from the data of INeill et al.l 
(2009) up to align it with our results. The trend due 
to p c is weaker than in observations, suggesting that p c 
contributes to the observed trend but that other effects 
also play a part. 

4. CONCLUSIONS AND FUTURE WORK 

We simulated a suite of 150 SNIa models with a range 
of p c to study the trends for a population of SNela. We 
find that o n average pro g enitor s with higher p c produce 
less 56 Ni. iHoflich et aU (|2010D argue that 56 Ni in the 
central regions of the exploding WD does not contribute 
to the light curve at maximum, and therefore they do 
not see a significant trend with central density in the 
maximum V-band magnitude, but rather in late-time 
brightness. The pure-deflagration models of lRopke et al.l 
( 2006) exhibit a shallow increase of produced 56 Ni as cen- 
tral density in c reases , in contradiction of our findings. 
Iwam oto et al.l (|1999h find that the trend with central 
density depends on the DDT transition density; extrap- 
olating from their results, our value of pddt should yield 
an increasing 56 Ni yield as central density increases. We 
find that small perturbations of the initial flame surface 
not only influence the final 56 Ni yield, but also its de- 
pendence on central density through variations in the 
duration of the deflagration phase caused by differences 
in plume development. The variation that follows from 
perturbations on the initial conditions is a critical aspect 
of multi-dimensional modeling. Only after many real- 
izations with different perturbations of the initial flame 
surface are simulated docs a statistically significant trend 
with central density emerge. This result, illustrated by 
Figure [21 demonstrates the need for an ensemble of sim- 
ulations to explore systematic effects in SNela. 

By relating p c to i coo i and 56 Ni to Ami5(B), our results 
support the observational finding that SNela from older 
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stellar populations are systematically dimmer. While a 
degeneracy between age and metallicity in the integrated 
light of stellar populations exists, the observed depen- 
dence of mean brightness of SNel a on mean stellar age 
is apparently the stronger effect (jGallagher et al.l 120081 : 
iHowell et al.l 12009) . Accordingly, our choice to neglect 
metallicity effects and consider only the effect of central 
density on 56 Ni yield allows us to offer a theoretical expla- 
nation for this observed trend. If we additionally consider 
the effect of metallicity, we may see a slightly stronger 
trend of decreasing brightn ess with increasing ag e as has 
been previously suggested (jTimmes et al.ll2003[ ). Other 
effects besides progenitor central density and metallicity, 
such as progenitor main sequence mass, may also con- 
tribute to this trend. 

The insensitivity of the overall Fe-group yield to cen- 
tral density, and therefore delay time, along with the 
dependence of the 56 Ni yield on central density, implies 
that SNela of similar brightnesses (and therefore simi- 
lar 56 Ni yield) from progenitors of different ages will not 
have the same total Fe-group yield. Those from older 
populations will, on average, have larger masses of stable 
species. This may argue for a slight non-u niformity in the 
Philli ps relation ba sed on environment ( Wooslc v et al.l 
I2007t IHoflich et al.|[2010l) . The resulting closely-related 
family of brightness-decline time relations also provides a 
physical motivation for intrinsic scatter in the Phillips re- 
lation as a result of combining populations with different 
mean stellar ages. In this picture the primary param- 
eter is the degree of expansion at DDT, determined by 
the morphology of the early flame (and the DDT den- 



sity, which we hold constant), and the age acts a weaker 
secondary parameter. In any case, the possibility of such 
an effect motivates further exploration of the impact of 
central density on the light curve itself. 
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